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We present a method to analyze magnetic properties of frustrated Ising spin models on 
specific hierarchical lattices with random dilution. Disorder is induced by dilution and ge- 
ometrical frustration rather than randomness in the internal couplings of the originated 
Hamiltonian. The two-dimensional model presented here possesses a macroscopic entropy 
at zero temperature in the large size limit, very close to the Pauling estimate for spin- ice 
on pyrochlore lattice, and a crossover towards a paramagnetic phase. The disorder due to 
dilution is taken into account by considering a replicated version of the recursion equations 
between partition functions at different lattice sizes. An analysis at first order in replica 
number allows for a systematic reorganization of the disorder configurations, leading to a re- 
currence scheme. This method is numerically implemented to evaluate the thermodynamical 
quantities such as specific heat and susceptibility in an external field. 

PACS numbers: 75.10.Hk, 05.50.+q, 75.10.Kt, 75.50.Lk 



I. INTRODUCTION 



Hierarchical lattices possess some interesting features of geometrical scale invariance, and sev- 
eral methods have been developed extensively in the past to study thermodynamical properties of 
spin models considered on such structures [l, 2|. The similarity with Migdal-Kadanoff renormal- 
ization procedure for magnetic spin systems 3] with exact recursion relations between coupling 
constants at the successive stages of lattice construction makes such models suitable for the study 
of critical properties. Indeed, phase transitions in these systems have non-mean field critical ex- 
ponents, effective dimensions and frustration effects within local loops. However, the absence of 
translation and the site dependent connectivity inherent to such theoretical construction make the 
physics apparently different from real systems. The importance of exact solutions in disordered 
lattices are highlighted by the possibility to obtain controllable recursion relations which gives 
interesting properties for Ising spin glass and quenched disorder models 0, Q|, or disordered Potts 



models U Mult icril ical point locat ions can moreover be cheeked carefully using Xishimori sym- 



metric lines 



9| and duality exact properties [101 ] . Frustration on hierarchical lattices has been 



studied for example on a diamond shape ferromagnetic structure with an additional transverse anti- 

n 

ferromagnetic coupling between two spins [ill ]. Important properties of the phase diagram when the 
coupling is increased concern the low temperature entropy which presents steps at specific values of 
this coupling reflecting the frustration character of the lattice geometry. These transitions are gen- 
erated by slight variations of the transverse coupling inducing a change in the nature of the ground 
states. Recursion relations can be found exactly in this non-disordered but frustrated model, hence 
equations for ground state structure, residual entropy and magnetization can be obtained using 
scaling properties of the hierarchical geometry. Experimental realizations of frustrated magnetic 
structures with spin dilution can be found in spin-ice materials doped with rare earth elements, such 
as Dy2- x Lua;Ti207, Dy2_ x Y :r Ti207, and Bx^-zYj^iO? [12 ]. where x parametrizes the random- 
ness [13j. These compounds are obtained from primary materials Dy2Ti2C>7 or Ho2Ti2C>7, where 
magnetic ions Dy 3+ and Ho 3+ occupy the sites of a pyrochlore lattice made of tetrahedra connected 
to each other by their vertices. The first magnetic rule for each individual tetrahedron is given by 
2 spins out of the tetrahedron and 2 spins in, known as first ice rule [ijj]. These magnetic ions are 
then replaced by non-magnetic atoms Y or Lu, giving rise to an experimental system of site dilution 
where a macroscopic fraction of sites is non magnetic. One of the interesting physical properties of 
these systems is the non-monotonic behavior of the zero-temperature entropy with respect with the 
dilution level 13] from the measured specific heat, which can be explained quite accurately within 



the Pauling approximation where tetrahedra are treated as independent 
with Monte-Carlo simulations on Husimi cactus 15[] is a step beyond the mode 



131 ] . Refined calculations 



of Pauling which 



U|. In the Bethe- 



takes into account the correlations between tetrahedra on the pyrochlore lattice 
Peierls approximation, a tree-like structure is constructed from a central site where the number of 
branches grows with the distance from the central site. However in this geometrical approach, the 
branches never reconnect together, a construction which forbids geometrical loops, but recursion 
equations can be in general written in a systematic manner. In this paper, we address the question 
of site disorder importance on thermal properties in the same spirit as 11( for pure systems. We 



indeed consider an antiferromagnetic tetrahedron-like and hierarchical structure with Ising spins 
cTj = ±1 located on the vertices with a probability < x < 1. Disorder configurations are given 
by a set of quenched random variables e« = 0, 1 located on each site i, as shown in Fig. Q] and 
starting from a single antiferromagnetic bond at level r = 0. The zero temperature entropy of the 
tetrahedron-like structure at level r = 1 can be evaluated exactly as we will see in the following and 
it appears to be non-monotonic as function of the dilution. We will then generalize the calculation 
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for larger sizes r > 1 using recursion equations in order to obtain information about frustration 
and disorder effects in the thermodynamical limit. This simple model can be considered as a good 
example on how to characterize spin-ice or spin-liquid states and how these states may emerge for 
example from the specific heat or susceptibility measurements. 

II. NOTATIONS AND METHOD 

We consider the hierarchical lattice displayed in Fig. Q] and constructed from a single antiferro- 
magnetic bond of Boltzmann strength K = J/T with coupling unity J = 1 at step r = 0. Each 
site on the bond can be vacant with probability < x < 1 or occupied by a spin particle with 
probability (1 — x). Tetrahedral like structure is replicated from the previous structure by adding 
at each step additional transverse antiferromagnetic bounds of value 5K > (see Fig. [T] (b) and 
(c)). In general, we use a set of variables e« = for site i, when the site is empty with probability 
i or £j = 1 when it is occupied with probability 1 — x, with Ising spins = ±1 on each occupied 
site. The total number of vertices N r at level r is equal to 

N = 2,N 1 = 4, N r+1 = 4N r -4,N r = ^ + a + 8 ■ (1) 

o 

Also the total number of bounds B r is equal to B r = 4 r + 2 r+1 — 2, and the maximum number of 
bonds connecting the two extreme sites is equal to L r = 2 r . We may define an effective dimension d 
such that B r = Lf, which gives d = 2 in the large size limit r ^> 1. In the following we will mainly 
concentrate on the frustrated case 5 = 1. Partition function for this dilute hierarchical lattice made 
of antiferromagnetic bonds can be constructed starting at step r = using the following standard 
method. Also an external magnetic field h e is taken into account in order to evaluate the spin 
susceptibility and magnetization. 

At recurrence r, we define the partial partition function (ea^e'a') which depends on some 
disorder configuration {r/} given by a set of vacancy distribution {e^ = 0, l}j with i corresponding 
to the sites located inside the structure delimited by the two extreme sites (ea, e'er'). These latter 
sites have a given disorder configuration {e, e'} which does not belong to {rf\. A recursion relation 
can be written between steps r and r + 1 of the lattice construction 

Z r { 5(ea,6V) = Tr^^^ex^Z^ 

x exp(- <5^ei<7ie 2 <7 2 + + e ^)) x expi-SKeae'a'^-^ 2 - 1 (2) 
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FIG. 1: Construction of the hierarchical lattice, starting from a single link K between two spins (a). 8K is 
the coupling between boundary sites, and the shaded areas in (c) define the previous lattice structures at 
level r — 1. 



where {rj} represents the disorder configuration {r]} = {rji, rj2, r)3, 774, e%, €2} and the Tqi sign 
the sum over all the possible spin values. The additional factors 2 £l ~ 1 2 e2 ~ 1 compensate the fact 
that we sum over eventual ghost spins located on vacant sites. The initial condition Fig. Q^a) is 
given by the two-spin partition function 



Z^ } {ea, eV) = exp(-lveeW) 



(3) 



which has actually no internal disorder dependence, {r/} = {0}. The magnetic field is imple- 
mented in Eq. ([2]) only for spins not belonging to the top and bottom vertices. Using replica 
method, we consider n copies of Eq. ([2]) for a given configuration {r/} and perform the average 



7 I J rr'< 



,a=l 



(4) 



Here [.]„ is meant for averaging over disorder {r/}. In this case, the recursion relation Eq. ([2]) 
becomes 



Z r+1 (ea a ,e'a a ') = exp(-5Kee' a a a' a ) J P(e 1 )P(e 2 )de 1 de 2 Jr {a? ^ 



} 



xZ r (ea a , eiaf )Z r (ei<rf , e'a' a )Z r (ea a , e 2 af)Z r (e 2 af , eV 1 



x exp 



rr/, , \ _ol _a 1 e \ V_a! , o\ r»n(ei — l)r,n(e2 — 1) 

<3Keie2 2^ cr i CJ 2 + ^ 2^( (J i + a 2) 2 2 



(5) 
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where P is the bimodal distribution for the dilute sites P(e) := x8(e) + (1 — x)5(e — 1). In the 
previous expression, we perform only the sum over the two spins af and connecting pairs of 
diamond shape structures. Starting with initial condition 

Z {ea a ,e'a a ') = exp(-Kee / ^ a a a' a ) (6) 

a 

all quantities Z r (ea a , e'a a ) can in principle be computed step by step. Free energy is 
evaluated directly after averaging over the remaining disorder and summation over the two spin 
variables. We define the complete partition function as 

z r (n) := I' P(e)P(e')dede'2^~ l h^'~ 1 hr Waya} Z r (ea a , eV*') exp + a' a )} (7) 

and take the limit n — > 



- KF {r) := lim - (z r (n) - 1) = 4(0) (8) 
n^O n 

after noticing that z r (0) = 1. The disorder averaged entropy can be expressed in terms of z' r (0) 
via the usual thermodynamical relations 

SW(T) := -— = <(0) - ^^4(0) (9) 

as well as the averaged specific heat = -Td 2 F^/dT 2 = K 2 d 2 z' r {0)/dK 2 . The relation 
between entropy and specific heat is given by the integral 

-YrdT' (io) 

which is a direct experimental way to measure the zero-temperature entropy by extrapolation, 
knowing that all spins are in the paramagnetic phase at high temperature, each contributing 
with a factor ln(2). The linear susceptibility corresponds to the excitation of the magnetic order 
parameter M = Yl l i a i under a field h e and is defined as = 9[< M >] v /dh e , after averaging 
over disorder. We can also relate to the free energy, using Tx^ = [< M 2 >] v — [< M > 2 } r , 
and [< M >} ri = —dF^ /dh e , in particular 

(r) d[< M >} v _Q^_ 

x dh e dh 2 • 1 ' 
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III. RESIDUAL ENTROPY FOR r = 1 STRUCTURE 

As a simple application, we consider the case r = 1 consisting in only 4 spins Fig. QJb), where all 
thermodynamical quantities can be computed exactly. The zero temperature entropy is identical 
to the entropy of a three-dimensional tetrahedron made of antiferromagnetic bonds. After iterating 
Eq. ([5]), we obtain 



Zi(e(rW Q ) = expt-Sifee'^VV 1 



a=l 



(1 - xfT Yl ( e SK + e~ SK cosh \lK{to a + e'a' a ) - 2 



-2z(l - x)2 n cosh 



K(ea a + eV Q ^ 



r 



(12) 



Then the function z\{n) is equal, after some computation, to 



Zl (n) = (l-x) 4 4™ 

+ 4x(l-x) 3 4 n 



1 + { cosh(4K) + e~ 4K sinh 2 (2^)}e- 2 ^ + e 2SK + cosh(2^) 
cosh(^)e 5K + { cosh(2i0 cosh(^) cosh(2^) 



sinh(2K) sinh(^) s inh(2^)}e 



-5K 



2r,n 



+ 2x 2 (l - x)*2 



+ 4x d (l -x)2 n cosh" (t^) + x 



cosh(2^)e-^ + e SK ) n + 2 ( cosh^e"* + e K " " 
h f 



(13) 



When K is large and in absence of external field, this expression can be put into the following 
form, by keeping the dominant contributions in the exponential terms 



n(s k -Ke k ) 



(14) 



This expansion is useful in order to identify the zero temperature entropy which can be written 
as a sum of contributions S^'(T = 0) = pkSk- Each partial entropy Sk and energy physically 
corresponds to a configuration of disorder with weight pk- In this limit, the asymptotic form for 
z\ (n) when < 5 < 1 is indeed equal to 



zi(n) ~ (l-x) 4 2V( 4M ) + 4x(l-x) 3 2V( 2 ^ 

+ 4x 2 (l - x) 2 2 n e n ^ + 2x 2 (l - x) 2 2 n e n5K + Ax 3 (l - x)2 n + x\ (15) 
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After taking the limit n = 0, the entropy and energy per spin are respectively equal to (T = 
0)/[4(l - x)\ = (1 - x 4 )/[4(l - a?)] ln(2) and E^ /[4(1 - x)] = -(1 - x) 3 (l - 5/2) - x{l - x) 2 {2 - 
5) - x 2 (l - x)(l + 6/2). The frustrated case 5 = 1 is particular since we obtain instead 



Zl (n) ~ {l-x) 4 6 n e n2K + Ax(l - x) 3 6 n e nK + 6x 2 {l - x) 2 2 n e nK 

+ 4x 3 (l -x)2 n + x A (16) 

and, in this specific case, the entropy and energy per spin have a different expression 

= -5(1 - - - "? - lAl - *)• (17) 

As expected, the entropy for a single tetrahedron in the frustrated case is larger. Plot of the 
entropy is given in Fig. [4] as function of x and appears to be non-monotonic with a local maximum 
around x = 0.295. The exact entropy value in absence of dilution is given by ln(6)/4 ~ 0.4479 
which is larger than for a system of coupled tetrahedra in the large size limit (see sections below). 

Susceptibility is derived exactly using definition Eq. (jlip and is plotted in Fig. [2] as function 
of the external field for different values of x. As expected, when x > 0, a peak appears at h e = 
when the probability to find a single spin in the tetrahedron is non zero. A simple calculation for 
the non-disordered case leads to three possible ground states per spin Eq depending on the value 
of the field: two spins up and two spins down with Eq = —2, three spins up and one spin down 
with Eq = —2h e , and all spins up with Eq = 6 — 4h e . Two transitions occur respectively at h e = 1 
and h e = 3. 

For larger sizes however, it is interesting to analyze the susceptibility using recursive equations 
Eq. ([5]) for the partition function Z r (a, a') in the absence of disorder x = 0. The recursive equations 
are detailed in appendix El The method is to assume an effective form for the partition function, 
with Ising effective coupling, magnetic field, and weight coefficient at all recursion levels. The 
exact recurrence equations for these three quantities as function of the external magnetic field h e 
allow us directly to compute the susceptibility. Plots of x^ r=6 \h e ) and magnetization per spin 
m ( r =^)(ji e ) are displayed in Fig. [3] showing multiple field transitions and magnetization plateaus 
up to h e = 127 for A% = 2732 spins. In the next section, the analysis is extended to the dilute 
case. We use the set of recursive equations Eq. ([5]) to compute the partition function at a given 
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FIG. 2: (Color online) Susceptibility per spin x^ 1 ' /[4(1 — x)] at T = 0.1 as function of the external field h e , 
for different values of dilution factor. In absence of disorder (black lines), the susceptibility is composed of 
two peaks located at h e = 1, 3. Inset: variation of the second peak position with dilution. 

step r. The main technique is to reorganize the partition function in the same manner as Eq. (|16p 



iterative procedure at low temperature on effective Ising constant couplings to obtain the limiting 
coupling distribution can in principle be performed, but not for the free energy, unlike the case of 
spin glass models. For example, for the disordered Potts model [a] where the disorder is given by a 
discrete distribution of couplings, zero temperature free energy satisfies some recursive equations. 
In general, when the disorder is located only in the couplings, the four structures represented as 
shaded areas in Fig. (He) have independent internal disorder, even if they are connected together 
by the four vertices. Therefore disorder can be factorized and direct iteration of Ising couplings by 
renormalization can be performed as well as the partition function weights such as I r in Eq. (|A1|) 
which are important for the determination of free energy and entropy. In the case of dilution 
however, the disorder configuration in the boundary vertex sites given by set {e, e' , e±, €2} is actually 
shared between the four diamond structures. For example variable ei is common to the two shaded 
diamonds on the left hand side of Fig. QJc). This makes the evaluation of weight coefficients I r 
Eq. (jAljl problematic in the disordered case since after few iterations correlations will develop 
rapidly. This is main reason why it is more convenient to consider a replica version of Eq. (|Aip and 



where 




One difficulty with dilution is that an 
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FIG. 3: (Color online) Susceptibility (blue line) X (^e) and magnetization m^ r \h e ) (red line) per spin 
at T = 0.05 as function of the external field h e in absence of disorder for level r = 6 (2732 sites). Inset: 
magnetization per spin in the large held limit until saturation. 



Eq. ([2]) with partial integration over the disorder located strictly inside the diamond structures, as 
explicitly defined previously by Eq. ([5|) . As we will see, an approximation scheme can be derived 
from the replica method and thermodynamical functions can be obtained. 



IV. RECURSION RELATIONS IN THE GENERAL CASE 

At step r, we assume from the previous analysis an expansion of the partition function in terms 
of configurational weights x k (l — x) nr ~ k C k r where n r is here the number of sites localized between 
the two extreme top and bottom sites and which satisfies the equation n r+ \ = 4n r + 2, with initial 
condition no = and N r = n r + 2. The main idea, as discussed before, is to keep a generic 
and minimal expression for the partial partition function, which depends on unknown functions 
satisfying recurrence equations such that 
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Z r {ea a ,e'a ,a ) = exp(-5Kee' ^ a a a' a ) >P x k {l - x) nr ~ k C k r exp (nl£\e, e 

k=0 



K^ee' £ A"> + ffM(e, e') + eV"]] 



(18) 



where H^\e, e') and i^(e, e') are symmetric functions of e and e'. Coupling if^ is independent 
of the boundary disorder. All these functions depend implicitly on if and /i e . As initial condition 
we impose if,^ = —if + 5if and iig°^(e, e') = Io°\e, e') = 0. We also take <5 = 1 in the following 
for the frustrated version. At the next level, Z r+ \ is evaluated by considering the product of four 
partition functions Z r as written in Eq. ([5]). 

Using the ansatz Eq. (fTBj) . this product is expanded as 



Z r+1 = exp(-5Kee'^2a a a' a ) ]T C£C^C^C*> fel+fe2+fe3+fe4 (l - ^-k.-^-k.-k, 

a ki,k2,k3,k4=0 

x J P(e l )P(e 2 )de l de 2 Jr {(T; ^ } exp (nl^(e, ex) + nljgfa, e') + nlg(e, e 2 ) + nl£(e 2 , e') 
x exp I -SKe^J^ a l a 2 - SK^i^i + ^)(e(T a + eW a )\ 2 n (^~ l h n ^-^ 

\ a a / 

a a « a / 

ex p (< E[- a + + < E[ eVa + + < } E^ Q + ^ + < + 

\ a a a a / 

[^5>i< + (19) 



x exp 



:h 

x exp 

The last term takes into account the missing field on former boundary spins erf and which 
are now summed up. It is useful to introduce the operator 1 = X^t=o ^fc,fci+fc 2 +fc3+A:4 or the integral 
form 

1 = E = E / r e^ ( - fe+fel+fa+fe+fc4) (20) 

k=0 k=0 J ° 71 

which is then inserted in the previous expression in order to reorganize the sum over the ki 
into a single sum over weights x k {l — x) inr ~ k C k nr . The integral over 6 can be performed using 
a first order expansion in n, sufficient to obtain 2^(0). For example, given integer p > 1 and 
< k < pn r , let us consider n r + 1 field variables ipi and define the quantity W n ^[(/Jj] made of 
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the product of p sums J2 k . C™£ exp(nyj^), i = 1, ■ • • ,p by analogy with the product of sums that 
appears in Eq. (|19|) . In addition, we impose the constraint Ya=i = ^ by using the Kronecker 
integral Eq. (|20p . and perform an expansion at first order in n 



W n> k[<Pl] ■ 



r 2n jo V ( n r 

JU i=l \ki=0 

r 2ir jo T> 

Jo 2vr 11 



fej e i0ki+ntp k . 



(21) 



2 *" d0 



l + np(l + e^)-"^C^e^V fcl 



n e -^(i + e i0 )^ 

o 27r 



min(fc,n r ) 

E 

fci=max(0,fc— (p— l)nr) 



nk\(~ik—k\ 

U nr U (p-l)n r 



Vfcl 



The exponentiation in the last line allows us, at first order in n, to reorganize the product of 
the p = 4 sums in Eq. ([5]) as a single sum over configurations of k vacant sites, with combinatorial 
factor Cp Ur . In particular, taking a constant value ipi := ip, we easily find W^fy] = Cp nr e nptp , and 
therefore the identity 



min(fc,n r ) 

E ' 

fcl=max(0,fc— (p— l)n r ) 



/~iki(~ik—k\ 
L ra T .°(p-l) nr 

nk 

pn r 



(22) 



from which we deduce that factors 



k,ki ._ ^ n r(p-l)n r 



n r ,p 



nk 

^pn r 



(23) 



can be considered as natural weights since the sum over integers k% is normalized. We may 
apply this technique to spin operators a a a' a appearing in Eq. (|19p as well, which are sum of n 
terms. Considering for example linear spin operator E^Q,o" a , and instead of Eq. (12 1 [> the function 



i=l \k t =0 

J ° ikA \i=l ) i 



(24) 
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We then take the limit n — > and define w^i] ■= dW n k[tpi]/dn\ n= o, so that 



'IT 

{ki} M=l 

Comparing this expression with the approximation 



2 cosh (if 



(25) 



™fcM - C pn r ln 



2cosh(p#^P£> fel 

ki 



(26) 



coming from the same analysis made in Eq. (|2ip . we can discuss two different cases. First, let 
choose <pi := cp, both Eq. (|25[) and Eq. (126p are identical, using normalization property Eq. (I22j) . 
Then, we may try non-constant fields such as ip\ := lip, which gives after summation the exact 
result 



(27) 

which is also identical to expression Eq. ([26]) using equality £V £>n^p^i = k/p. Considering 
Eq. (|2ip for spin operators will be useful to simplify expression Eq. (|19p and obtain recursive 
equations in presence of dilution. Computing recursive equations for (e, e'), 1^ (e, e'), and ii^ 
follows two steps: the integral over variable 9 coming from the constraint k\ + ki + ks + = k, 
with k{ = 0, • • • ,n r , introduced in Eq. (|19p . is performed using transformation Eq. (|21|) over the 
different spin operators with p = 4. This will allow for the partial summation over spins af and > 
and averaging over random variables e\ and e2, see Fig.QJc). We obtain an expression for Z r+ \ as a 
summation over configurations x k+l (l — x ) Anr+2 ~ k ~' ] 'C\ nr C\, with 4n r sites coming from inside the 
shaded diamond structures in Fig. HHc), plus the two sites coming from the integration over ei and 
€2- This sum can be furthermore reorganized using again identity Eq. (12ip in order to finally obtain 
Eq. (fl~8P at level r + 1 as a sum over weights x k (l — x) inr+2 ~ k C k nr+2 (with 4n r + 2 = n r+ i) and 
new coupling values. All details of this development are presented in appendix [Bl where recurrence 
equations are written explicitly in Eq. (IBlip . 

Using ansatz Eq. (|18p . and integrating over the boundary site degrees of freedom, the complete 
partition function Eq. (|7|) is then equal to 



w k [(pi = lip] 



2 cos\\{kKip) 
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■(») = £c*x fc (l-x)"' 



:i-x) 



fc=0 



exp(n4 r - ) (l,l))2 w {e< -^cosh [2^(1,1) + 2^ 



+ 2s(l - x) exp(n4 r) (0, l))2 n cosh n [#£ r) (0, 1) + ^] + x 1 exp(n4 r) (0, 0)) 
and the free energy is derived directly from the previous equation 



+ e" 



-K 



KF& = z'M = Y J Ctx k {l-xr- k [{l-xf{H2) + lt ) (1,1) 



+ In O 



fc=0 

^' r) - 5x cosh 



2#£ r) (l,l)+2^ 



+ e k 



K ( „ r) +SK 



+ 2x(l - x) { ln(2) + j£° (0, 1) + In cosh ) (0, 1) + ^] } + x 2 /^ (0, 0)] . (28) 

After a few steps, the number of configurations is growing rapidly, as the number of weights 
C k r x k (l — x) nr becomes exponentially large as well as the number of iterative functions to 
evaluate. We can obtain however a very good approximation is we notice that these weights are 
distributed closely around a Gaussian when n r is sufficiently large 



(fc— xn r — x+i) 2 

C k n X k (\ - x]^- k ~ ^ \ 2x(l-x >^_ (29) 

\J 2irx(l — x)n r 

For r = 4 for example, the number of internal sites is equal to n r = 170, and the previous 
approximation is very accurate. Numerically, we solved the iterative functions up to level r = 4 
included, using Eq. (jBlip . and then apply for higher levels r > 4 the Gaussian approximation for 
k distributed with 4 standard deviations around the mean value xn r + x — 1/2, which gives precise 
results. 



V. CALORIMETRY AND THERMODYNAMICAL FUNCTIONS IN THE DILUTE 

CASE 

In this section, we evaluate different thermodynamical quantities as function of dilution using 
Eq. (|28p . The residual entropy per spin is plotted as function of x in Fig. [J] for r between 1 and 
6. For the single tetrahedron structure r = 1, the entropy is numerically identical to the exact 
expression Eq. (|17|) and presents non-monotonic dependence with increasing dilution. For r larger, 
the entropy is reduced, but saturates rapidly after r = 5 which corresponds to 684 sites. In the limit 
of extreme dilution, the entropy per spin is simply equal to ln(2) as expected. It is interesting to 
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FIG. 4: (Color Online) Residual entropy per spin S^/[N r (l-x)] at zero temperature computed from the 
expression of the free energy Eq. (|28l) . The entropy limit per spin for a very dilute system (independent 
spins) is close to ln(2) = 0.6931 as expected. The orange dashed line represents the Pauling residual entropy 
Sp(x) given in Eq. (f30|) . with a value for the undilute case equal to Sp(x = 0) = \ ln(3/2) ~ 0.2027 (see 
text). Inset: zero temperature entropy for the undilute system as function of the inverse system size. First 
exact values are = ln(6)/4 = 0.4479, = ln(2)/3 = 0.2310, S^/U ~ 0.2297. The term at 

recursion level 6 is approximately equal to S'( 6 )/2732 ~ 0.20804. The dashed straight line is the Pauling 
entropy Sp(0). 

compare the resulting entropy with the Pauling estimation Sp for an infinite number of tetrahedra 
treated as independent as function of dilution [13] 



S P (x) = ln(2) - 3x 2 (l - x) ln(2) - 2x(l - x) 2 ln(4/3) - ^(1 - x) 3 ln(8/3). (30) 

Comparison between Sp and shows very similar values at low dilution, especially the entropy 
difference for the undilute case is quite small, 5 P (0) = |ln(3/2) ~ 0.2027 and S^/N 6 ^ 0.20804 
which is an upper bound (see also inset of Fig. 0]). Exact values for the undilute hierarchical 
structure can be computed up to a certain order but the entropy shows a behavior similar to spin 
ice models. Approximations on pyrochlore lattice made of Ising antiferromagnet tetrahedra give a 
closer value around 0.20410 11711 . 
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FIG. 5: (Color online) Specific heat per spin C v r /[N r (l — x)] at zero field and level r — 6 (2732 sites) for 
different disorder probability values. Specific heat for the pure case x = is derived from the exact recursion 
equations Eq. (|A4|) given in appendix [A] The dashed blue line at x — 0.95 is the fit with a two-level model 
which accounts for the Schottky anomaly at a temperature close to unity, see text and Eq. pip . Inset: low 
temperature behavior where a plateau is visible at x = 0.3. The dashed black line indicates the position of 
each Schottky peak as function of dilution. 

(r) i—i 

The specific heat Cy is displayed in Fig.[5jas function of temperature for five different values of 
x. The curves presents in general a broad maximum or Schottky anomaly at a temperature around 
T = 1 corresponding to the typical coupling J = 1 and associated with a crossover between a low 
temperature spin-ice state and paramagnetic state. The system however stays antiferromagnetic 
in the low temperature regime but is highly degenerated. The ground state energy per spin can be 
computed exactly for the first terms in absence of dilution /4 = —1/2, E&> /12 = —5/6, and 
the limiting value is estimated to be I?( r>1 ) /N r ~ —0.9 using the recurrence equations in appendix 
lAl At large dilution, the specific heat can accurately be fitted with a two-level model with gap A 
and constant Co, as it can be seen in Fig. [5] 



A 2 e A / T 

c a PP rox = c (31) 

For example, the curve for x = 0.95 was fitted with the previous formula using Cq ~ 0.114 
and A = 1.988, which corresponds to the specific heat for a gas of dilute pairs of spins with 
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FIG. 6: (Color online) Low temperature behavior of specific heat per spin C { v r) /[N r (l-x)T 2 } at zero field 
and level r — 6 (2732 sites) for different disorder probability values. A local maximum is developing for low 
dilution 0.1 < x < 0.3, then non exponential behavior is observed for intermediate values x ~ 0.5. 

Co — 2(1 — x) = 0.1 and energy coupling very close to J = A/2 = 1. At lower temperature 
however T ~ 0.1, the specific heat presents a second broad peak at intermediate dilution factor 
x ~ 0.3 (see inset of Fig. [5] and Fig. [6|) which can not be reproduced by a two-level model. 

These characteristics are exemplified in Fig. [6] for T < 0.1. The exponential-like Arrhenius 
behavior of the pure system seems to evolve to more complex features associated with a very small 
peak contribution at intermediate dilutions and non-exponential deviations. Arrhenius behavior is 
then recovered when we approach large dilution modeled by Eq. (|3ip . We have actually rescaled 
the specific heat in Fig. [6] by a factor 1/T 2 , in order to check if a phonon-like two-dimensional 

(r) rl 

Debye contribution is seen at intermediate dilution, where Cv should scale like T with d = 2 in 
our model and where /T 2 should saturate in this scaling regime. Such scaling was analyzed for 
example in pyrochlore compound Bi2Ti2C>7 (with d = 3) in order to measure the excess of specific 
heat due to Einstein oscillator additional contributions that could give rise to a broad peak at low 
temperature [l8|]. The scaling in T 2 in Fig. [6] is more appropriate since a T 3 scaling would present 
clearly a divergence. To summarize, the specific heat per spin as function of both temperature 
T and dilution factor x is displayed in Fig. [TJ where the variation of the main Schottky peak 
amplitude with x shows a decreasing behavior towards a system made of individual pairs of spins 
with a broader extension. 




FIG. 7: (Color online) Surface plot of the specific heat per spin Ci, r> /[N r (l-x)] at zero field and level r = 6 
(2732 sites) as function of temperature and dilution factor x. The Schottky peak amplitude is reduced as x 
increases. 

Susceptibility curves as function of dilution and field are plotted in Fig. [8j We chose to represent 
only the low-field excitations h e < 1 in order to follow the displacement and amplitude of the first 
peaks with dilution, in particular those corresponding to Fig. [3J in the same low field region. As 
dilution is increased, a new peak appears at h e = corresponding to excitations of uncoupled and 
isolated spins. The location of the peak at h e = 1/2 does not change except its amplitude. The 
peak located at h e = 1/3 is instead moving towards higher field values, with several intermediate 
peaks in the range 1/5 < h e < 1/2. Smaller peaks at h e < 1/3 are moving towards the origin 
instead. A surface plot Fig. Ogives a general view of how peaks are moving with field with respect 
of dilution, and how their amplitude vanishes as we approach the high dilution regime. 
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FIG. 8: (Color online) Behavior of the susceptibility x /[N r (l — x)] at low temperature (T=5.10 -2 ) for 
2732 sites (r — 6) as function of the field h e and for several dilution factors x. Here are represented only 
low field excitations h e < 0.65. The values of x are plotted on the left axis with a color corresponding to 
each curve for clarity. 

VI. CONCLUSION 



In this paper we propose a method to study frustrated hierarchical lattices in presence of 
dilution based on the replica method and reorganization of configurational weights at first order 
in the replica parameter n. Interesting properties of the dilute spin- ice state at low temperature 
can be examined within this approximation by implementing recursive equations for the partition 
function and leading to specific heat and susceptibility as function of temperature and external 
field. Clear crossover evidence is seen between spin-ice and paramagnetic states in the specific 
heat with the presence of a Schottky peak, and the zero-temperature entropy follows closely the 
Pauling approximation at least at moderate dilution. Specific heat presents also a secondary 
peak at low dilution (x ~ 0.3) and very low temperature with non-Arrhenius behavior. This 
feature is probability due to a combination of low Einstein and phonon-like modes. This makes 
this hierarchical model a good candidate for exploring in details the physics of spin-ices or spin- 
liquids. Additional analysis can probably be made using correlation functions for example or short 



range ferromagnetic order parameter 



111 ] to probe the spin correlations in the low temperature 



state. This approximation scheme based on replica may probably be implemented more easily to 



19 



n 4 
3.5 




O.I 0.2 0.3 0.4 0.5 0.6 

h e 
(b) 

FIG. 9: (Color online) (a) Surface plot representing the spin susceptibility x /[-^r(l — %)] at low temperature 
(T=5.10 -2 ) for 2732 sites (r — 6) as function of the field h e and dilution factors < x < 1. Here are 
represented only low field excitations 0.1 < h e < 0.65 corresponding to the first three peaks of Fig. [3] (b) 
Map view of the surface plot. 

hierarchical spin glass models with modal distribution of couplings, since the quenched disorder 
is treated as independent between the recursive diamond structures, making the need of a partial 
partition function not necessary and therefore simplifying the analytical recurrence. We would like 
to acknowledge M. Gingras for useful discussions on thermal properties in spin-ice systems. 
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Appendix A: Recursion relations for the non-disordered model 

In this appendix we write the recursion relations for the non-disordered case {x = 0). Starting 
from the partition function Zq(ct, a') = exp(-Kaa') of a single antiferromagnetic link between two 
spins a and a', and 5 = 1, we can generally assume the following recursive and stable form at any 
step r 



Z r {a, a') = exp [l r — K r aa' + H r (a + a') 



(Al) 



with Iq = 0, Kq = K, and Hq = as initial conditions. At the next level r + 1, we form the 
product 



Z r+ i(a,a) = Jr^ aicr , 2 jZ r (a,ai)Z r (a 1 ,a')Z r (a,a2)Z r (a2,cr') 
x exp ^ - Kaa - Kaia 2 + -^r [°i + °2]) • 



(A2) 



After replacing Z r by its ansatz Eq. (|Alj) . and performing the sum over the internal spin degrees 
of freedom, we obtain 



Z r+ i{a,a') = exp (4I r - Kaa' + 2i/ r (a + </) ) <^ 2cosh 2K r (a + <r') - 4ff r - 2 



e-^ + 2e^ 



The term into bracket {• • • } can be rewritten as 



2K r (a + </) - 4F r - 2-j; e~ K + 2e^ = exp ( I r - K r oo + H r {a + a'; 



2 cosh 

with the following equations for I r , K r and H, 



(A3) 



exp(/ r — + 2H r ) = 2 cosh 
exp(I r — K r — 2H r ) = 2 cosh 
exp(/ r + K r ) = 2 cosh 



AK r - AH r - 2H r 



4K r + 4H r + 2H, 
4K r + 2H 



e- K + 2e K , 



e~ K + 2e K . 



This set of equations can be solved by eliminating successively the arguments in the exponential 
terms. Then the recursive solutions for the new couplings of Z r+ \ are given by 
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'r+l 



4J r + I r , K r+ i = K + K r , H r+ i = 2H r + H r . (A4) 



These relations can be easily implemented in order to compute numerically the different ther- 
modynamical quantities from free energy F r = —TlnZ r . 

Appendix B: Recursion relations at finite temperature 



In this section the recursive equations for the different couplings in Eq. (I18j) are derived. Using 
Kronecker integral Eq. (|20p . we can rewrite Z r+ \ as 



Z r+l = Y^x k (l-xf^- k f P(e l )P(e 2 )de l de 2 Jr { ^^ ] 2 n ^- l h<^ (Bl) 

fc=0 

x exp ( -5Kee' a a a' a - 5K ei e 2 ^ a? a% - 5K ^(eiof + e 2 a%)(ea a + e'a' a ) ] 

\ a a a / 

r h 1 f^W rift Ur 

x ^ [y E^< + ^ )] / Z e ~ t9k E ^.C^C n fc 3C^e ie (^ + ^ +fc 3+^) 

o fci,/c2,fc 3 ,A:4=0 

x exp {ni£>( Cl £1 ) + ^(61,0 + ^(e, 6 2 ) + n^(6 2! e') 

+ <W E ^ + V E + E ^ + *£V E 

a. a a a 

+ tf Cl ) J> CT a + £l of ) + < } ( e ', £1 ) £>V Q + £l o?) 

a ct 

+ fl£>( e , e 2 ) ^><r Q + e 2 a%) + < } (e' , e 2 ) £>V Q + e 2 crf ) j . 



Using Eq. (|2T|) . we can integrate over 9 and rewrite Eq. (|B1|) as 
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4n r „ 

Y,CL r x k (l-x)^- k / P(e 1 )P(6 2 )de 1 de 2 Tr KjCT . } 2"^- 1 )2™^- 1 ) 

fc=0 

exp ( -5Kee' ^ a Q cj ,Q - <Sifeie 2 E <r? - <5if E( eiCT i + e 2 <7f )(e<7 Q + eV a ) 

\ a a a 



min(?i r ,fc) 

fci=max(0,fc— 3n r ) 
min(ri r ,fe) 

+ E ' <4< ) E(^ a + e '0(^?+^) 

fci=max(0,fc— 3n r ) ol 
min(n r ,fc) 

+ E ^; fc 4E[< } ( e ' e i)( OTa + e i<) + < } ( e/ ' e i)( eVa + e i<) 

fci=max(0,A;— 3n r ) <* 

+ pM( e , e2 )( e(J a + e2 a 2 Q ) + < ) (e',6 2 )(6V a + 6 2 a 2 a )l + ^ £( Cl o? + ea o£) 



(B2) 



Now the sum over intermediate spins erf and o"2 can be performed directly. Let first define 
intermediate couplings 



K 



(r) 



min(ra r ,fc) 

E <l< ] - S*> 

fcl=max(0,fc— 3n r ) 



(B3) 



and new functions 



min(n r ,fc) 

fci=max(0,fc— 3n r ) 
min(n r ,fe) 

fci=max(0,fc— 3n r ) 

Then we isolate the part in Eq. ()B2p containing only <rf and (Xg , and perform the sum: 



(B4) 



+ 



Tr {(Tf >B o } exp f -ffifeie 2 E + £< r) E( e i< + e 2 CT 2 )(^° + eV a ) 

\ « a 

[Jff« (e, e>, ei) + ^]E £ i< + ( e > ^ '*) + ^] E 

a a / 

\{ [exp(-^ ei e 2 )2co S h [K^(e 1 + e 2 )(ea a + eV Q ) + ^\e, e', e x ) + ^)e x + (i^ r) (e, e', e 2 ) + 



exp(^ ei e 2 )2cosh ^ r) ( £l - e 2 )(ea a + eV to ) + (#f (e,e' )ei ) + ^)e 1 - (tff(e,e',e 2 ) + ^Je 2 



23 



Next, we perform the integration over e\ and €2 



Z r+1 = J]cy(l-i)^exp -<51fee'J><V to 

k=0 V a / 

(1 - x) 2 exp (nl k (e, e', 1, 1) + i^ r) (l, 1, 1) ^(ea a + eV to )^ 
x 2"[] |exp(-«5K) cosh [2K^ (ea a + e'a' a ) + 2^ r) (e, e', 1) + 2^] + exp(5K) 

+ 2x(l-x)exp ^n4(e,e',0,l) + ^ r) (l,0,l)^(ea a + eV Q )^ 



x 2 n J] cosh [i^ r) (e^ + eV Q ) + i^ r) (e, 6', 1) + ^ 
+ x 2 exp ( nl k (e, e', 0, 0) + #£ r) (0,0,1)^ (ea a + eV 



(B6) 



We introduce now a set of functions {M^^e, e'), Qk,h H k ,i{e, e')} for each of the three terms 
appearing in the previous expression and proportional to (1 — x) 2 (I = 0), 2x(l — x) (/ = 1), and 
x 2 (I = 2) respectively. The first factor associated with (1 — x) 2 can be exponentiated such that 



exp (4(6, e', 1, 1) + ff^fl, 1, l)(ea a + eV to )) 



x 2 <J exp(-<5.FQ cosh 2K ( k r) (ea a + e'a' a ) + 2H^ , (e, e' , 1) + 2^ +exp(^) 



(r). 



T 



= : exp (M M (e, 6') + Q M eeW a + H kfi (e, e')(ea a + eV a )) . 
Similarly we have for the second term proportional to 2x(l — x) 



(B7) 



exp (l fc (e,e',0, 1) + i^ r) (l,0, l)(ea a + e'a' a )) 2 cosh \k { k r] \ea a + eV to ) + H { k r \e,e\ 1) 



T . 

= : exp (M M (e, e') + Q M eeVV a + fT M (e, e')(ea a + e'a' a )) . (B8) 

The last term associated with x 2 can be exponentiated using the values M^^e, e') := 
4(e,e',0,0), Q fci2 := and H K2 {e,e') := H { k \o,0, 1). All these functions can be identified in 
a unique way by using the four possible configurations for a a and a' a . Then the partition function 
can be rewritten as 
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4n r 



Z r+1 = Cl r x k (l - x)^~ k exp -SKee' £ a a a' a (B9) 

;=0 V a / 

' 2 

£ C\x\\ - xf~ l exp (nM kjl (e, e') + Q^ee' £ a a a' a + if*,, (e, e') ^(e<r a + e'a' a )) 

.1=0 ol a 

As before we can expand the exponential terms at first order in n and rearrange the powers in 
x such that 



(\ 4n r +2=n r+ i 
-5Kee' £ aV a £ C 4 fc nr+2 x fc (l - x) 4 "-+ 2 ^ x (BIO) 
a / fc=0 



ex P ( E [nM fe _ M (e, e') + Q k _ hl ee> £ + H k _ l:l (e, e') ^(ea a + eV to )] ) . 

Z=max(0,fc-4n r ) 4n r +2 Q a 

From this result, we can deduce finally the recursive equation for the couplings 



min(2,fc) r k-lril 

A .<,. +1 , _ ^ ^1Q„. V , 

Z=max(0,fc-4n r ) 4n r +2 
min(2,fe) ^k-ls-il 

lt +l \^')= E ^M fc _ M ( e ,e'), and 

Z=max(0,fc-4re r ) °4n r +2 



min(2,fc) pk-lpl 
Z=max(0,fc-4n r ) 4n r +2 
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